% This file is used to define the Global variables, initial conditions and call the required files.
% The variables are:
% p: represents zeta : wave height
% q : represents phi : surface velocity potential
% 2-D case .. in X direction


%--- main parameters
%clear all;

global Ux Uy  h  g  p  q r tmax Lx Ly dx dy Nx Ny y tmax;
%----- Physical Parameters----------

Ux = 1; % Current Velocity in the X direction. Assumed Constant
Uy=  1;  % Current velocity in the Y direction
h = 50; % height of the bed, assumed constant
g=  9.81; % Gravitational Constant


%------------ Simulation Parameters % -----------------
tmax = 1;     % max time, s
Lx=  100;         % Length in the X direction
Ly =100;            % Length in the Y direction
Nx =20;         % number of axial nodes in the X direction
Ny=20;          % Number of Nodes in the Y direction
steps_x = Nx-1;    % number of axial steps
steps_y = Ny-1;
dx = Lx/steps_x;   % axial step spacing, delta x
dy = Ly/steps_y;

% Check the CFL condition:
dt1 = 0.5*dx/sqrt(g*h);
%--- Initial Conditions ----------
pi1=1;
qi=1;
r=zeros(Nx*Ny,1);
for j=1:Nx
    for i=1:Ny
        index = (j-1)*Ny+i;
        p(index,1) = pi1*cos(4*pi*(j-1)*dx/Lx);
        q(index,1) = qi;
        r(index,1)=  1.0;
    end
end


%Run the matlab benchmark Code
i=1;
Benchmark(2);
A(1:Nx*Ny,1)=y(end,1:Nx*Ny);
B(1:Nx*Ny,1)=y(end,Nx*Ny+1:2*Nx*Ny);
R(1:Nx*Ny,1)=y(end,2*Nx*Ny+1:3*Nx*Ny);


%dt_v =[0.1 0.01 ] ;
%for i=1:length(dt_v)
  %  dt=dt_v(i);
 %   TimeDiscretizers(dt,1);
 %   A(1:Nx*Ny,i+1)=p;
 %   B(1:Nx*Ny,i+1)=q;
   % C(1:Nx*Ny,i+1)=r;
  %  X = sprintf('The Norm of the Error for Zeta is %f' , norm_L2(A(:,1),A(:,i+1)));
  %  disp(X);
  %  X = sprintf('The Norm of the Error for varphi is %f' ,  norm_L2(B(:,1),B(:,i+1)));
  %  disp(X);
%end

%for k=1:length(dt_v)+1
k=1;
    for j=1:Nx
        for i=1:Ny
            index = (j-1)*Ny+i;
            p_surf(i,j,k) = A(index,k);
        end
    end

%end
figure;
hold on
colormap([1 0 0;0 0 1])
[xx yy] =meshgrid(1:Nx,1:Ny);
surf(xx,yy,p_surf(:,:),ones(Nx));
surf(xx,yy,m_Zeta(:,:),ones(Nx)+1);
%surf(xx,yy,m_Z1(:,:),ones(Nx)+2);
%surf(xx,yy,p_surf(:,:,1),ones(Nx));
%surf(xx,yy,p_surf(:,:,3),ones(Nx)+1);
view(17,22)
% 
% % 
% 
% for j=1:Nx
%     for i=1:Ny
%         index = (j-1)*Ny+i;
%         p_surf(i,j) = p(index);
%     end
% end

